function dataset = generate_counterfactual_mto(Prob, Ss, rs, lambdas, tol, max_It, rs_series, thresh);
    
global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha...
        xi xi1 ub lb pi_e H w theta Int_pi_w...
        sigmae elasts Pr log_e_mu log_e_sigma skill_shock  phi sigma gamma rts fp r_bar r_base w_mult shares...
        pe thetas tmp_mult adj w_bar ss_calc unemp h_resid shock_multiplier myslope eta mu_a thetas myiter part_eqm  thresh_flag
    
    thresh_flag = thresh; % Imports the scale of the MTO policy (income threshold)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial Choices Before MTO
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Pre-allocate and compute variables
    b0 =  b ;
    N_n = length(Ss);    
    eta = 1;

    % Empty matrix for data
    dataset = cell(5, 4);
    ln_e_hat = -1/2*log_e_sigma^2; % mean of the non-normalized idiosyncratic wage shock
    
    % Calculate choice probabilities
    [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
    
    % Initial Distribution
    myDist0 = P_S .* Prob;
    myDist0  = sum(bsxfun(@times, myDist0 , reshape(pi_theta, [1,1,1,2])),4);
    myDist1 = myDist0;
    
    N_size = squeeze(sum(myDist0,1:2));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate statistics pre MTO
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    t = 1;
    tax = 0;

    % Income pdf and cdf
    pi_w = squeeze(sum(Prob, 2:3));
    cdf_w = cumsum(pi_w,1);
 
    % Gini
    my_prod = pi_w .* pi_w' .* abs(w - w');
    gini = sum(my_prod(:))/(2*dot(pi_w, w));
    rich_cutoff = [0.8];
    poor = rich_cutoff(1); 
    rich = 1 - rich_cutoff(1); 
    
    % Average education, ability, wage
    avg_ed = squeeze(sum(educ_mat.* myDist1,1:2)./sum( myDist1 ,1:2))';
    avg_ability = squeeze(sum(a'.* myDist1,1:2)./sum( myDist1 ,1:2))';
    avg_wage    = squeeze(sum(sum(w .* myDist1,1),2)./sum(sum( myDist1 ,1),2))';

    % Dissimilarity
    P_poor_S = zeros(length(Ss),1);
    P_rich_S = zeros(length(Ss),1);
    
    find_rich_w = @(x)interp1(w,cdf_w,x,[],'extrap')-0.8;
    w_rich = fzero(find_rich_w,(w_max+w_min)/2);
    agg_cdf = cumsum(squeeze(sum(myDist1, 2)), 1);
        
    for j_c = 1:length(Ss)
        P_poor_S(j_c) = P_poor_S(j_c) + interp1(w,agg_cdf(:, j_c)',w_rich);
        t1 = interp1(w,agg_cdf(:,j_c)',max(w)) - interp1(w,agg_cdf(:,j_c)',w_rich);
        P_rich_S(j_c) = P_rich_S(j_c) + t1;
    end
    dissim = 0.5 * sum(abs(P_poor_S ./ poor - P_rich_S ./ rich)) ;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Rank-Rank Correlation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% Upward Mobility
    % find 25th percentile (cdf_w was already defined above)
    N_bp1 = 125;
    q = linspace(0,1,N_bp1);
    w_q = zeros(1,N_bp1);
    w_q(1) = w_min;
    w_q(N_bp1)= w_max+1e-5;
    for i = 1:(N_bp1-2)
      find_q_w = @(x)interp1(w,cdf_w,x,[],'extrap')-q(i+1);
      w_q(i+1) = fzero(find_q_w,(w_max+w_min)/2);
    end

    w_qk = zeros(1,N_bp1);
    w_qk(1) = w_min;
    w_qk(N_bp1)= w_max+1e-5;
    cdf_w_kids = cumsum(sum(Prob,2));
    for i = 1:(N_bp1-2)
      find_q_w = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-q(i+1);
      w_qk(i+1) = fzero(find_q_w,(w_max+w_min)/2);
    end

    parents_rank = zeros(N_w,N_a);
    for m = 1:(N_bp1-1)
        for i = 1:N_w
            if w_q(m)<=w(i) && w_q(m+1)>w(i)
                parents_rank(i,:) = m;
            end
        end
    end

    wprime = unemp + (b0 + F_es) .* f_w(w);
    wprime = repmat(wprime, [1 1 1 N_e 2]);
    tmp = repmat(e, [1 N_w N_a N_n 2]);
    tmp = permute(tmp, [2 3 4 1 5]);

    wprime = wprime .* tmp;

    [~, kids_rank] = histc(wprime,w_qk);
    pi_joint_w_a_e = repmat(myDist1,1,1,1,N_e) .* permute(repmat(pi_e, [1 N_w N_a N_n]), [2 3 4 1]);
    mu_prank = sum(sum(sum(sum(parents_rank.*pi_joint_w_a_e))));
    mu_krank = sum(sum(sum(sum(kids_rank.*pi_joint_w_a_e))));
    mu_krank = pi_theta*mu_krank(:);
    prank_std = sum(sum(sum(sum(sum((parents_rank-mu_prank).^2.*pi_joint_w_a_e)))));
    prank_std = sqrt(prank_std(:));
    krank_std = sum(sum(sum(sum(sum((kids_rank-mu_krank).^2.*pi_joint_w_a_e)))));
    krank_std = sqrt(pi_theta*krank_std(:));
    rank_cov = sum(sum(sum(sum(sum((parents_rank - mu_prank).*(kids_rank-mu_krank).*pi_joint_w_a_e)))));
    rank_cov = pi_theta*rank_cov(:);

    rank_corr = rank_cov/(prank_std*krank_std);

    % Q1 to Q1
    [~,idx] = min(abs(cdf_w -.25));
    find_25_w = @(x)interp1(w,cdf_w,x,[],'extrap')-0.25;

    pi_w_kids = sum(Prob,2:3);
    cdf_w_kids = cumsum(pi_w_kids);
    find_25_w_kids = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-0.25;
    w_25_kids = fzero(find_25_w_kids,(w_max+w_min)/2);

    % find probability of staying below 25th percentile for each (ability,wage) pair
    Q1toQ = zeros(N_a,N_w, N_n);

    for t2=1:2
        for i=1:N_a
        for j=1:idx
                for c = 1:N_n
                    my_share = Prob(j,i)*P_S(j,i, c,t2)*pi_theta(t2);
                    wj = w(j);
                        Q1toQ(i,j) = Q1toQ(i,j) + ...
                            my_share*sum(pi_e((b0 + F_es(j,i,c)).*e.*f_w(wj) <= w_25_kids));
                end
            end
        end
    end
    mobility_up = sum(Q1toQ(:))/cdf_w(idx);
    

    dataset{t,1} = dissim;
    dataset{t,2} = gini;
    dataset{t,3} = Ss';
    dataset{t,4} = N_size';
    dataset{t,5} = rs';
    dataset{t,6} = r_base;
    dataset{t,7} = avg_ability;
    dataset{t,8} = rank_corr;
    dataset{t,9} = avg_ed;
    dataset{t,10} = P_poor_S' ./squeeze(sum( myDist1 ,1:2))';
    dataset{t,11} = tax;            
    dataset{t,12} = mobility_up;
    dataset{t,13} = 0;
    dataset{t,14} = 0;

    N_size = squeeze(sum(myDist0,1:2));
    denom  = sum(myDist0 .* wprimeS, 1:3);
    bet1   = bet1 / (denom)^xi1;
    Ss     = Ss .* (denom)^xi1;
    
    % Welfare
    Prob_1  = P_S .* Prob;
    Prob_1  = sum(bsxfun(@times, Prob_1 , reshape(pi_theta, [1,1,1,2])),4);
    [P_S_1, ~, wprimeS, ~] = choiceProb_ct_par(Ss, rs, 0, 0, Ss);
    [~, ~, ~, ~, Eu]       = choiceProb_ct_par_eu(Ss, rs, Prob_1, 0, 0, Ss);
    myDist_1  = sum(bsxfun(@times, P_S_1 .* Prob_1 , reshape(pi_theta, [1,1,1,1,2])),5);
      
    % Movers statistics
    avg_wage_c = sum(myDist1 .* wprimeS, 1:5);
    percent_gain = 0;
    move_mass = 0;
    for c = 1:N_n
        counter_neigh = setdiff(1:length(Ss), c);
        mymass = sum(myDist_1(:,:,c,counter_neigh) ,4);
        move_mass = move_mass + squeeze(sum(mymass,1:3));
        percent_gain = percent_gain + sum(mymass .* (wprimeS(:,:,c,:)), 1:3);
    end
    percent_gain = squeeze(percent_gain) ./ move_mass;

    c = 3;
    move_mass = 0;
    percent_gain = 0;
    mymass = sum(myDist_1(:,:,c,1) ,4);
    move_mass = move_mass + squeeze(sum(mymass,1:3));
    percent_gain = percent_gain + sum(mymass .* (log(wprimeS(:,:,c,1)) - log(wprimeS(:,:,c,c))), 1:3);
    percent_gain_ng = squeeze(percent_gain) ./ move_mass;
    
    move_mass = 0;
    percent_gain = 0;
    mymass = sum(myDist_1(:,:,c,1) ,4);
    move_mass = move_mass + squeeze(sum(mymass,1:3));
    percent_gain = percent_gain + sum(mymass .* (log(wprimeS(:,:,c,1)) - log(wprimeS(:,:,c,c))), 1:3);
    percent_gain_geo = squeeze(percent_gain) ./ move_mass;

    dataset{t,15} = squeeze(sum(myDist_1,1:2)) ;
    dataset{t,16} = Eu ;
    dataset{t,17} = percent_gain_ng;
    dataset{t,18} = percent_gain_geo;
    dataset{t,31} = squeeze(sum(sum(a'.*myDist1,1),2))./squeeze(sum(sum(myDist1,1),2)) ;
    dataset{t,32} = avg_wage;
    
    % Invert the housing equation to get lambdas
    lambdas =  inv_hs(shares , rs); 
    r_base = dot(w, sum(Prob,2:3));
    
    % Update the distribution to the following period
    w_max0 = max(w);
    w_min0 = min(w);
    
   Prob = P_S .* Prob; % Neighborhood of parent's choice
   Prob  = sum(bsxfun(@times, Prob , reshape(pi_theta, [1,1,1,2])),4);
   tax = 0; w_mto = 0;
   myiter = t;
   Probm = update_dist_int_edu_mto(Prob, Pr, Ss, rs, tax, w_mto, lambdas, 1);
   F_es = skill_shock .* F_es; % Adjust parent's wage for shock
   Prob = Probm;

   [P_S, ~, ~, ~] = choiceProb_ct_par(Ss, rs, 0, 0, Ss);
   myDist0 = P_S .* Prob;
   myDist0  = sum(bsxfun(@times, myDist0 , reshape(pi_theta, [1,1,1,1,2])),5);
    
   % Load steady state spillovers and rents for partial equilibrium
   % exercises

    load Ss_ss
    load rs_ss
    load tax_ss
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Dynamics> Introduce MTO experiment
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for t = 1:5
        myiter = t+1;

        if t==1
            % Reshape distribution
            myDist1 = P_S .* Prob;
            myDist1  = sum(bsxfun(@times, myDist1 , reshape(pi_theta, [1,1,1,1,2])),5);
            dataset{t,15} = squeeze(sum(sum(myDist1,1),2)) ;
            dataset{t,31} = squeeze(sum(sum(a'.*myDist1,1),2))./squeeze(sum(sum(myDist1,1),2)) ;
            id_max=[1 2 3];
        end

        
        
        % Generate shock to spillover
        if t == 1
            eta = skill_shock;
        end
        if t >= 1
            pi_w = squeeze(sum(Prob, 2:3)); % Wage PDF
            cdf_w = cumsum(pi_w); % Wage CDF
            % Calculate MTO threshold
            [~, idx] = min(abs(cumsum(pi_w) - thresh));
            w_mto = w(idx);
            idx_mto = idx;

            % Calculate thresholds for control groups
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.0015));
            w_ct  = w(idx_ct);
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.05));
            w_ct5  = w(idx_ct);
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.10));
            w_ct10  = w(idx_ct);
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.15));
            w_ct15  = w(idx_ct);
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.20));
            w_ct20  = w(idx_ct);
            [~, idx_ct] = min(abs(cumsum(pi_w) - 0.25));
            w_ct25  = w(idx_ct);
        end
        
        r_base = dot(w, sum(Prob,2:3));
        w_bar = r_base;

        w0 = w;
        w_max0 = w_max;
        w_min0 = w_min;
        thresh_grid = [0.0015 0.05 0.10 0.15 0.20 0.25];
                    
         if part_eqm == 1
            % In partial equilibirum, use steady state spillovers and rents
            Ss=Ss_ss(:,t+1);
            rs=rs_ss(:,t+1);
         else
             % In general equilibrium, calculate spillovers and rents
             [Ss, rs] = find_spillovers_mto(Ss, lambdas, rs, Prob, tol, max_It, w_mto);
         end
   
        
        % Calculate tax and mean wage
        [~, tax, mean_wage] = find_rents_mto(rs, Ss, lambdas, Prob, w_mto);
        
        % In partial equilibrium set the tax to 0
        if part_eqm == 1
            
            tax = tax_ss(find(thresh_grid==thresh),t+1);
        end
        
        % Identify neighborhoods by mean wage
        [~, id_max] = sort(-mean_wage);
        
        % Neighborhood Choice
        [P_S, F_es, wprimeS, educ_mat] = choiceProb_ct_par(Ss, rs, tax, w_mto, mean_wage);
        [~, ~, ~, ~, Eu,Eu_n,Eu_mat,Eu_i] = choiceProb_ct_par_eu(Ss, rs, Prob, tax, w_mto, mean_wage);
        pdf = Prob;
        choices = P_S;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Calculate Statistics
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        N_size = squeeze(sum(bsxfun(@times, Prob .* P_S , reshape(pi_theta, [1,1,1,1,2])),[1,2,3,5]));
       
        myDist1 = P_S .* Prob;
        myDist1  = sum(bsxfun(@times, myDist1 , reshape(pi_theta, [1,1,1,1,2])),5);
        
        % Share of individuals eligible for the voucher
        vouch_share = sum(sum(sum(sum(myDist1(:,:,id_max(end),:).* (w <= w_mto)) )));
        
        % Share of individuals who take up the voucher
        take_up = sum(sum(myDist1(:,:,id_max(end),id_max(1)).* (w < w_mto)) )/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_mto)) ));
       
        % Flag individuals who take up the voucher
        take_up_id = myDist1(:,:,id_max(end),id_max(1)).* (w < w_mto);
        take_up_id2=take_up_id;
        take_up_id(take_up_id>0) = 1;
        take_up_id(take_up_id==0) = 0;
        
        % Expected Kids' Income for MTO and control group
        exp_inc       = myDist1.* wprimeS;
        income_mto    = sum(sum(exp_inc(:,:,id_max(end),id_max(1)).* (w < w_mto)) )/sum(sum(myDist1(:,:,id_max(end),id_max(1)).* (w < w_mto)) );
        income_mto_ct = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct)) ));
        income_mto_ct5 = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct5)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct5)) ));
        income_mto_ct10 = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct10)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct10)) ));
        income_mto_ct15 = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct15)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct15)) ));
        income_mto_ct20 = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct20)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct20)) ));
        income_mto_ct25 = sum(sum(sum(exp_inc(:,:,id_max(end),:).* (w < w_ct25)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct25)) ));
        
        income_pol    = sum(sum(sum(exp_inc(:,:,id_max(end),:)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:)) ));

        N_s_ct   = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct)) ))';
        N_s_ct5  = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct5)) ))';
        N_s_ct10  = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct10)) ))';
        N_s_ct15 = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct15)) ))';
        N_s_ct20 = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct20)) ))';
        N_s_ct25 = squeeze(sum(sum(myDist1(:,:,1,:).* (w < w_ct25)) ))';

        
        % Education for MTO and control group
        educ = educ_mat.* myDist1;
        educ_mto = sum(sum(educ(:,:,id_max(end),id_max(1)).* (w < w_mto)) )/sum(sum(myDist1(:,:,id_max(end),id_max(1)).* (w < w_mto)) );
        educ_mto_ct   = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct)) ));    
        educ_mto_ct5  = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct5)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct5)))) ;
        educ_mto_ct10 = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct10)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct10))) );
        educ_mto_ct15 = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct15))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct15)))) ;
        educ_mto_ct20 = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct20))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct20)))) ;
        educ_mto_ct25 = sum(sum(sum(educ(:,:,id_max(end),:).* (w < w_ct25)) ))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct25))));


        ab      = a'.* myDist1;
        ab_ct   = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct)) ));    
        ab_ct5  = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct5))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct5)))) ;
        ab_ct10 = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct10))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct10))) );
        ab_ct15 = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct15))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct15)))) ;
        ab_ct20 = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct20))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct20)))) ;
        ab_ct25 = sum(sum(sum(ab(:,:,id_max(end),:).* (w < w_ct25))))/sum(sum(sum(myDist1(:,:,id_max(end),:).* (w < w_ct25))));
        
        curr_ab = myDist1.*a';
       
        curr_ability_mto   = sum(sum(curr_ab(:,:,id_max(end),id_max(1)).* (w < w_mto)) )/sum(sum(myDist1(:,:,id_max(end),id_max(1)).* (w < w_mto)) );
        curr_ability_nomto = sum(sum(sum(curr_ab(:,:,id_max(end),id_max(2):id_max(3)).* (w < w_mto))) )/sum(sum(sum(myDist1(:,:,id_max(end),id_max(2):id_max(3)).* (w < w_mto)) ));
        
        %%% Gini %%%

        pi_w = squeeze(sum(Prob, 2:3));
        cdf_w = cumsum(pi_w,1);
       
        my_prod = pi_w .* pi_w' .* abs(w0 - w0');
        gini = sum(my_prod(:))/(2*dot(pi_w, w0)); 
       
        
        %%% Dissimilarity %%%
        rich_cutoff = [0.8];
        poor = rich_cutoff(1); 
        rich = 1 - rich_cutoff(1); 
        
        agg_cdf = cumsum(myDist1, 1);
        P_poor_S = zeros(length(Ss),1);
        P_rich_S = zeros(length(Ss),1);

        find_rich_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-0.8;
        w_rich = fzero(find_rich_w,(w_max+w_min)/2);
        for j_s = 1:N_n
            for j_a = 1:N_a
                for j_c = 1:length(Ss)
                    P_poor_S(j_c) = P_poor_S(j_c) + interp1(w0,agg_cdf(:,j_a,j_s, j_c)',w_rich);
                    t1 = interp1(w0,agg_cdf(:,j_a,j_s,j_c)',max(w0)) - interp1(w0,agg_cdf(:,j_a,j_s,j_c)',w_rich);
                    P_rich_S(j_c) = P_rich_S(j_c) + t1;
                end
            end
        end
        dissim = 0.5 * sum(abs(P_poor_S ./ poor - P_rich_S ./ rich)) ;

        % Average ability, education, wage
        avg_ability = squeeze(sum(a'.* myDist1,1:3)./sum( myDist1 ,1:3))';
        avg_ed      = squeeze(sum(educ_mat.* myDist1,1:3)./sum( myDist1 ,1:3))';
        avg_wage    = squeeze(sum(sum(sum(w.*myDist1,1),2),3)./sum(sum(sum( myDist1 ,1),2),3))';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        

        dataset{t+1,1} = dissim;
        dataset{t+1,2} = gini;
        dataset{t+1,3} = Ss';
        dataset{t+1,4} = N_size';
        dataset{t+1,5} = rs';
        dataset{t+1,6} = r_base;
        dataset{t+1,7} = avg_ability;
         dataset{t+1,9} = avg_ed;
        dataset{t+1,10} = P_poor_S' ./squeeze(sum( myDist1 ,1:3))';
        dataset{t+1,11} = tax;
        dataset{t+1,13} = vouch_share;
        dataset{t+1,14} = take_up ;
        dataset{t+1,15} = squeeze(sum(myDist1,1:2)) ;
        dataset{t+1,16} = Eu ;
        dataset{t+1,19} = income_mto;
        dataset{t+1,20} = [income_mto_ct income_mto_ct5 income_mto_ct10 income_mto_ct15 income_mto_ct20 income_mto_ct25];
        dataset{t+1,29} = [curr_ability_mto curr_ability_nomto];
        dataset{t+1,31} = squeeze(sum(sum(a'.*myDist1,1),2))./squeeze(sum(sum(myDist1,1),2)) ;
        dataset{t+1,32} = avg_wage;
        dataset{t+1,33} = Eu_n;
        dataset{t+1,34} = educ_mto;
        dataset{t+1,35} = [educ_mto_ct educ_mto_ct5 educ_mto_ct10 educ_mto_ct15 educ_mto_ct20 educ_mto_ct25];
        dataset{t+1,36} = [N_s_ct N_s_ct5 N_s_ct10 N_s_ct15 N_s_ct20 N_s_ct25];
        dataset{t+1,37} = [ab_ct ab_ct5 ab_ct10 ab_ct15 ab_ct20 ab_ct25];
        dataset{t+1,38} = income_pol;
        dataset{t+1,40} = take_up_id;
        dataset{t+1,41} = Eu_mat;
        dataset{t+1,42} = Eu_i;
        dataset{t+1,43} = pdf;
        dataset{t+1,44} = choices;
        dataset{t+1,45} = idx_mto;
        dataset{t+1,46} = wprimeS;
        dataset{t+1,49} = educ_mat;





        Prob = Probm;
        myDist0 = myDist1;
        
    end

end